Supplementary Figure 1. Hex plot of MD age effects

  • load age effect files and set up dataframes
scalar = "dti_md"
HCPD_ageeffects <- read.csv(paste0(output_root, "/", "HCPD", "/GAM/", scalar, "/HCPD_GAM_dev_measures.csv"))
HBN_ageeffects <- read.csv(paste0(output_root, "/", "HBN", "/GAM/", scalar, "/HBN_GAM_dev_measures.csv"))
PNC_ageeffects <- read.csv(paste0(output_root, "/", "PNC", "/GAM/", scalar, "/PNC_GAM_dev_measures.csv"))
# make a list of the different df names: [dataset]_GAM_ageeffects_[scalar]
datasets <- c("HCPD", "HBN", "PNC")
ageeffect_df_names <- c()
for (dataset in datasets) {
  ageeffect_df_names <- append(ageeffect_df_names, paste0(dataset, "_ageeffects"))
}

# format age effect df's
ageeffect_dfs <- lapply(ageeffect_df_names, format_ageeffect)
ageeffect.fdr_dfs <- lapply(ageeffect_dfs, sig_nodes)
## There are 2392/2400 significant nodes
## There are 2350/2400 significant nodes
## There are 2273/2400 significant nodes
names(ageeffect.fdr_dfs) <- ageeffect_df_names
tracts <- setdiff(unique(ageeffect.fdr_dfs$HCPD_ageeffects$tract_label), "Corticospinal") # CST for supplement
cor_HCPD_HBN <- round(cor.test(ageeffect.fdr_dfs$HCPD_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$HBN_ageeffects$GAM.smooth.AdjRsq)$estimate, 2)
cor_HCPD_PNC <- round(cor.test(ageeffect.fdr_dfs$HCPD_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$PNC_ageeffects$GAM.smooth.AdjRsq)$estimate, 2)
cor_HBN_PNC <- round(cor.test(ageeffect.fdr_dfs$HBN_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$PNC_ageeffects$GAM.smooth.AdjRsq)$estimate, 2)

# permute for p-values
HCPD_HBN_ageeffect_perm <- perm_cor_test(ageeffect.fdr_dfs$HCPD_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$HBN_ageeffects$GAM.smooth.AdjRsq)
HCPD_PNC_ageeffect_perm <- perm_cor_test(ageeffect.fdr_dfs$HCPD_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$PNC_ageeffects$GAM.smooth.AdjRsq)
PNC_HBN_ageeffect_perm <- perm_cor_test(ageeffect.fdr_dfs$PNC_ageeffects$GAM.smooth.AdjRsq, ageeffect.fdr_dfs$HBN_ageeffects$GAM.smooth.AdjRsq)

HCPD_to_corr <- ageeffect.fdr_dfs$HCPD_ageeffects %>% mutate(HCPD = abs(GAM.smooth.AdjRsq)) %>% mutate(tract_node = paste0(tractID, "_", nodeID)) %>% select(tract_node, HCPD)
HBN_to_corr <- ageeffect.fdr_dfs$HBN_ageeffects %>% mutate(HBN = abs(GAM.smooth.AdjRsq)) %>% mutate(tract_node = paste0(tractID, "_", nodeID)) %>% select(tract_node, HBN)
PNC_to_corr <- ageeffect.fdr_dfs$PNC_ageeffects %>% mutate(PNC = abs(GAM.smooth.AdjRsq)) %>% mutate(tract_node = paste0(tractID, "_", nodeID)) %>% select(tract_node, PNC)

HCPD_HBN_ageeffects <- merge(HCPD_to_corr, HBN_to_corr, by = "tract_node")
HCPD_PNC_ageeffects <- merge(HCPD_to_corr, PNC_to_corr, by = "tract_node")
HBN_PNC_ageeffects <- merge(HBN_to_corr, PNC_to_corr, by = "tract_node")

text_HCPD_HBN <- bquote(italic(r) == .(cor_HCPD_HBN) * "," ~ italic(p[perm]) < .(format(round(HCPD_HBN_ageeffect_perm$p_value, 4), scientific = FALSE)))
HCPD_HBN_ageeffects_plot <- hex_plot(HCPD_HBN_ageeffects, "HCPD", "HBN", text_HCPD_HBN, ylim1 = 0 , ylim2 = 0.64, xlim1 = 0, xlim2 = 0.62, x_text = 0, y_text = 0.64)

text_HCPD_PNC <- bquote(italic(r) == .(cor_HCPD_PNC) * "," ~ italic(p[perm]) < .(format(round(HCPD_PNC_ageeffect_perm$p_value, 4), scientific = FALSE)))
HCPD_PNC_ageeffects_plot <- hex_plot(HCPD_PNC_ageeffects, "HCPD", "PNC", text_HCPD_PNC, ylim1 = 0 , ylim2 = 0.64, xlim1 = 0, xlim2 = 0.62, x_text = 0, y_text = 0.64)

text_HBN_PNC <- bquote(italic(r) == .(cor_HBN_PNC) * "," ~ italic(p[perm]) < .(format(round(PNC_HBN_ageeffect_perm$p_value, 4), scientific = FALSE)))
HBN_PNC_ageeffects_plot <- hex_plot(HBN_PNC_ageeffects, "HBN", "PNC", text_HBN_PNC, ylim1 = 0 , ylim2 = 0.64, xlim1 = 0, xlim2 = 0.62, x_text = 0, y_text = 0.64)

# make bin colors the same across plots
# get the maximum bin counts from all datasets
max_count_HCPD_HBN <- max(ggplot_build(HCPD_HBN_ageeffects_plot)$data[[1]]$count, na.rm = TRUE)
max_count_HCPD_PNC <- max(ggplot_build(HCPD_PNC_ageeffects_plot)$data[[1]]$count, na.rm = TRUE)
max_count_HBN_PNC  <- max(ggplot_build(HBN_PNC_ageeffects_plot)$data[[1]]$count, na.rm = TRUE)
global_max_count <- max(max_count_HCPD_HBN, max_count_HCPD_PNC, max_count_HBN_PNC)
scale_limits <- c(0, global_max_count)
HCPD_HBN_ageeffects_plot <- HCPD_HBN_ageeffects_plot +
  paletteer::scale_fill_paletteer_c("ggthemes::Blue-Green Sequential", 
                                    direction = -1, 
                                    limits = scale_limits, 
                                    oob = scales::squish)
HCPD_PNC_ageeffects_plot <- HCPD_PNC_ageeffects_plot +
  paletteer::scale_fill_paletteer_c("ggthemes::Blue-Green Sequential", 
                                    direction = -1, 
                                    limits = scale_limits, 
                                    oob = scales::squish)

HBN_PNC_ageeffects_plot <- HBN_PNC_ageeffects_plot +
  paletteer::scale_fill_paletteer_c("ggthemes::Blue-Green Sequential", 
                                    direction = -1, 
                                    limits = scale_limits, 
                                    oob = scales::squish)


ageeffects_plots <- ggarrange(HCPD_PNC_ageeffects_plot, HBN_PNC_ageeffects_plot, HCPD_HBN_ageeffects_plot, nrow = 1, common.legend=T, legend = c("bottom"), labels = "auto")

title.grob <- textGrob("Age Effect Magnitude of Mean Diffusivity Across Datasets", 
                   gp=gpar(col="black", fontsize = 24, face = "bold"), hjust = 0.5)
 
grid.arrange(arrangeGrob(ageeffects_plots, top = title.grob))

#ggsave(paste0(suppfig_dir, "supp_hexplots.png"), grid.arrange(arrangeGrob(ageeffects_plots, top = title.grob)), height = 7, width = 20, units = "in")

Supplementary Table 1. Coefficient of Variation for tract ends and deep tract regions

  • load tract profiles
PNC <- readRDS(sprintf("%1$s/%2$s/tract_profiles/tract_profiles_for_NEST.RData", output_root, "PNC"))
HCPD <- readRDS(sprintf("%1$s/%2$s/tract_profiles/tract_profiles_for_NEST.RData", output_root, "HCPD"))
HBN <- readRDS(sprintf("%1$s/%2$s/tract_profiles/tract_profiles_for_NEST.RData", output_root, "HBN"))

tract profiles coefficient of variance

datasets <- list(HCPD = HCPD, HBN = HBN, PNC = PNC)

# create summary dfs
for (dataset_name in names(datasets)) {
  dataset <- datasets[[dataset_name]]
  make_summary_dfs("dti_md", dataset_name, dataset) # this creates summary_[dataset]_[scalar] dataframes that include mean, sd, se  
}
clipEnds = 5
bin_size = 5

PNC_cv <- summary_PNC_dti_md %>% 
  filter((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) | 
           (nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size)) | 
           (nodeID > (50-bin_size) & nodeID <= (50+bin_size))) %>%
  mutate(node_position = case_when((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) |
                                     nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size) ~ "superficial",
                                   (nodeID > (50-bin_size) & nodeID <= (50+bin_size)) ~ "deep")) %>%
  group_by(tract_label, tractID, node_position) %>%
  summarise(average_cv = round(mean(cv),2)) %>% 
  pivot_wider(names_from = node_position, values_from = average_cv, names_prefix = "average_cv_") %>%
  select(tractID, average_cv_deep, average_cv_superficial) %>% 
  rename(`PNC Deep` = average_cv_deep,
         `PNC Superficial` = average_cv_superficial,
         Tract = tractID)  
 
PNC_cv$Tract <- gsub("_", " ", PNC_cv$Tract)
PNC_cv <- PNC_cv[,-1]
PNC_cv <- PNC_cv %>% ungroup()

HCPD_cv <- summary_HCPD_dti_md %>% 
  filter((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) | 
           (nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size)) | 
           (nodeID > (50-bin_size) & nodeID <= (50+bin_size))) %>%
  mutate(node_position = case_when((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) |
                                     nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size) ~ "superficial",
                                   (nodeID > (50-bin_size) & nodeID <= (50+bin_size)) ~ "deep")) %>%
  group_by(tract_label, tractID, node_position) %>%
  summarise(average_cv = round(mean(cv),2)) %>% 
  pivot_wider(names_from = node_position, values_from = average_cv, names_prefix = "average_cv_") %>%
  select(tractID, average_cv_deep, average_cv_superficial) %>% 
  rename(`HCPD Deep` = average_cv_deep,
         `HCPD Superficial` = average_cv_superficial,
         Tract = tractID)  
 
HCPD_cv$Tract <- gsub("_", " ", HCPD_cv$Tract)
HCPD_cv <- HCPD_cv[,-1]
HCPD_cv <- HCPD_cv %>% ungroup()

HBN_cv <- summary_HBN_dti_md %>% 
  filter((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) | 
           (nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size)) | 
           (nodeID > (50-bin_size) & nodeID <= (50+bin_size))) %>%
  mutate(node_position = case_when((nodeID < (clipEnds + bin_size) & nodeID >= (bin_size)) |
                                     nodeID <= (99-clipEnds) & nodeID > (99-clipEnds-bin_size) ~ "superficial",
                                   (nodeID > (50-bin_size) & nodeID <= (50+bin_size)) ~ "deep")) %>%
  group_by(tract_label, tractID, node_position) %>%
  summarise(average_cv = round(mean(cv),2)) %>% 
  pivot_wider(names_from = node_position, values_from = average_cv, names_prefix = "average_cv_") %>%
  select(tractID, average_cv_deep, average_cv_superficial) %>% 
  rename(`HBN Deep` = average_cv_deep,
         `HBN Superficial` = average_cv_superficial,
         Tract = tractID)  
 
HBN_cv$Tract <- gsub("_", " ", HBN_cv$Tract)
HBN_cv <- HBN_cv[,-1]
HBN_cv <- HBN_cv %>% ungroup()

combined_cv <- left_join(PNC_cv, HCPD_cv)
combined_cv <- left_join(combined_cv, HBN_cv)
colnames(combined_cv) <- sub("^[^ ]+ ", "", colnames(combined_cv))
combined_cv$Tract <- gsub("Left", "L", combined_cv$Tract)
combined_cv$Tract <- gsub("Right", "R", combined_cv$Tract)

combined_cv <- rbind(combined_cv, c("Average Across Tracts", round(colMeans(combined_cv[,c(2:7)]), 2)))


cv_table_html <- combined_cv %>%
  kbl(format = "html", align = "lcccccc", digits = 2) %>%
  add_header_above(c(" " = 1, "PNC" = 2, "HCP-D" = 2, "HBN" = 2)) %>%
  add_header_above(c(" " = 1, "Average Coefficient of Variation" = 6), bold = TRUE) %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(0:nrow(combined_cv), extra_css = "line-height: 0.5;")  %>%
  column_spec(c(2, 4, 6), background = "gray90", include_thead = TRUE)
 
cv_table_html
Average Coefficient of Variation
PNC
HCP-D
HBN
Tract Deep Superficial Deep Superficial Deep Superficial
L Arcuate 4.26 4.11 3.8 3.87 7.83 5.25
R Arcuate 4.99 4.15 4.15 3.79 7.26 5.48
Callosum Anterior Frontal 5.06 3.53 4.41 3.54 8.29 4.83
Callosum Motor 11.42 4.63 8.34 3.77 10.97 7.04
Callosum Occipital 7.36 4.34 4.43 4.25 14.79 5.7
Callosum Orbital 6.98 4.08 4.23 3.87 12.04 6.37
Callosum Posterior Parietal 5.56 4.87 3.9 3.74 12.97 6.28
Callosum Superior Frontal 10.38 4.76 6.69 3.84 8.8 6.25
Callosum Superior Parietal 8.13 4.23 4.52 3.48 11.65 6.32
Callosum Temporal 8.61 4.76 5.53 3.82 14.68 9.45
L Corticospinal 3.1 10.84 3.22 6.24 7.33 8.74
R Corticospinal 3.4 10.73 2.87 6.56 7.02 8.86
L Inferior Fronto.occipital 9.67 4.58 6.01 4.2 7.59 5.48
R Inferior Fronto.occipital 5.55 4.95 5.32 4.37 7.16 5.49
L Inferior Longitudinal 4.79 4.66 4.09 3.93 6.5 6.43
R Inferior Longitudinal 5.7 4.57 4.14 3.81 6.87 5.62
L Posterior Arcuate 4.97 4.21 4.22 3.66 6.1 4.37
R Posterior Arcuate 5.19 4.59 4.1 3.65 6.22 4.9
L Superior Longitudinal 4.25 4.35 3.78 3.89 7.56 5.22
R Superior Longitudinal 4.49 4.35 3.46 3.81 7.35 5.43
L Uncinate 3.45 4.51 4.38 4.13 7.2 6.91
R Uncinate 3.29 4.27 3.44 4.14 6.62 6.37
L Vertical Occipital 4.77 3.9 3.88 4.26 5.48 5.05
R Vertical Occipital 5.65 4.03 3.94 4.05 5.57 5.18
Average Across Tracts 5.88 4.92 4.45 4.11 8.49 6.13 
#save_kable(cv_table_html, paste0(suppfig_dir, "table1_cv.html"))
#webshot(paste0(suppfig_dir, "table1_cv.html"), paste0(suppfig_dir, "table1_cv.pdf"))
#webshot(paste0(suppfig_dir, "table1_cv.html"), paste0(suppfig_dir, "table1_cv.png"))

Supplementary Figure 2. Age effect plots for corticospinal tract

color_HCPD = "#3FB8AFFF"
color_HBN = "#FF9E9DFF"
color_PNC = "#cc0468"
cst_PNC <- readRDS("/cbica/projects/luo_wm_dev/two_axes_manuscript/output/PNC/NEST/dti_md/Corticospinal_bin5_clip5_1sided.RData")
NEST_PNC <- data.frame(c("Corticospinal", cst_PNC$pval, "PNC")) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))

cst_HCPD <- readRDS("/cbica/projects/luo_wm_dev/two_axes_manuscript/output/HCPD/NEST/dti_md/Corticospinal_bin5_clip5_1sided.RData")
NEST_HCPD <- data.frame(c("Corticospinal", cst_HCPD$pval, "HCP-D")) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))

cst_HBN <- readRDS("/cbica/projects/luo_wm_dev/two_axes_manuscript/output/HBN/NEST/dti_md/Corticospinal_bin5_clip5_1sided.RData")
NEST_HBN <- data.frame(c("Corticospinal", cst_HBN$pval, "HBN")) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))
ageeffect_measure = "GAM.smooth.AdjRsq"
colors <- c("#3FB8AFFF", "#FF9E9DFF", "#cc0468")
names(colors) <- c("HCP-D", "HBN", "PNC")
bin_num_nodes = 5
clipEnds = 5

# prepare the data for making lollipop plots
df_all <- bind_rows(
  ageeffect.fdr_dfs$HCPD_ageeffects %>% mutate(Dataset = "HCP-D"),
  ageeffect.fdr_dfs$HBN_ageeffects  %>% mutate(Dataset = "HBN"),
  ageeffect.fdr_dfs$PNC_ageeffects  %>% mutate(Dataset = "PNC")
)

df_all$Dataset <- factor(df_all$Dataset, levels = c("HCP-D", "PNC", "HBN"))
df_formatted <- format_deep_superficial(df_all, ageeffect_measure, bin_num_nodes, clipEnds)  
df_formatted <- df_formatted[complete.cases(df_formatted), ]
df_formatted <- df_formatted %>% filter(tract_label == "Corticospinal")

# merge with NEST results
NEST <- rbind(NEST_HCPD, NEST_HBN, NEST_PNC)
lollipop_data <- prepare_lollipop_data(df_formatted)   
lollipop_data$Dataset <- factor(lollipop_data$Dataset, levels = c("PNC", "HCP-D", "HBN"))

# plot lollipops
lollipop_plot_CST <- make_lollipop_plot(tract = "Corticospinal", data = lollipop_data)
ageeffect_measure = "GAM.smooth.AdjRsq"
ageeffects_plot_CST <- plot_age_effect(tract = "Corticospinal", scalar = "dti_md", ageeffect_measure = ageeffect_measure, 
                                color_HCPD = color_HCPD, color_HBN = color_HBN, color_PNC = color_PNC, 
                           clipEnds = 5, ylim2 = 0.41, ylim1 = 0, fontsize = 18, legend_position = "none")
CST_plot <- plot_grid(ggdraw() + draw_label("Corticospinal", size = 18, y = 0.2, hjust = 0.5),
                                  plot_grid(plotlist = list(ageeffects_plot_CST, lollipop_plot_CST),
                                            align = "h", ncol = 2, rel_widths = c(1, 0.4)), ncol = 1, rel_heights = c(0.25, 1)) +  bgcolor("white") 

CST_plot

#ggsave(paste0(suppfig_dir, "supp_CST_ageeffect.svg"), CST_plot, height = 5, width = 5, units = "in")
#ggsave(paste0(suppfig_dir, "supp_CST_ageeffect.png"), CST_plot, height = 5, width = 5, units = "in")

Supplementary Figure 3. Age effect plots for fractional anisotropy (FA).

scalar = "dti_fa"
PNC_ageeffects <- read.csv(paste0(output_root, "/", "PNC", "/GAM/", scalar, "/PNC_GAM_dev_measures.csv"))
HCPD_ageeffects <- read.csv(paste0(output_root, "/", "HCPD", "/GAM/", scalar, "/HCPD_GAM_dev_measures.csv"))
HBN_ageeffects <- read.csv(paste0(output_root, "/", "HBN", "/GAM/", scalar, "/HBN_GAM_dev_measures.csv"))
# make a list of the different df names: [dataset]_GAM_ageeffects_[scalar]
datasets <- c("PNC", "HCPD", "HBN")
ageeffect_df_names <- c()
for (dataset in datasets) {
  ageeffect_df_names <- append(ageeffect_df_names, paste0(dataset, "_ageeffects"))
}

# format age effect df's
ageeffect_dfs <- lapply(ageeffect_df_names, format_ageeffect)
ageeffect.fdr_dfs <- lapply(ageeffect_dfs, sig_nodes)
## There are 1621/2400 significant nodes
## There are 1061/2400 significant nodes
## There are 1947/2400 significant nodes
names(ageeffect.fdr_dfs) <- ageeffect_df_names 
tracts <- setdiff(unique(ageeffect.fdr_dfs$HCPD_ageeffects$tract_label), "Corticospinal") # CST for supplement

# PNC: There are 1621/2400 significant nodes - 67.5%
# HCP-D: There are 1061/2400 significant nodes - 44.2%
# HBN: There are 1947/2400 significant nodes - 81.1%
ageeffect_measure = "GAM.smooth.AdjRsq"
ageeffects_plots <- lapply(tracts, plot_age_effect, scalar = "dti_fa", ageeffect_measure = ageeffect_measure, 
                                color_HCPD = color_HCPD, color_HBN = color_HBN, color_PNC = color_PNC, clipEnds = 5, ylim2 = 0.25, ylim1 = 0, fontsize
                           = 20, legend_position = "none")
names(ageeffects_plots) <- tracts

arrange callosal plots

arrange_callosum_plots_nolollipop(ageeffects_plots) + bgcolor("white")

#ggsave(paste0(suppfig_dir, "supp_ageeffects_fig1_FA.svg"), arrange_callosum_plots_nolollipop(ageeffects_plots), height = 10, width = 20, units = "in")
#ggsave(paste0(suppfig_dir, "supp_ageeffects_fig1_FA.png"), arrange_callosum_plots_nolollipop(ageeffects_plots), height = 10, width = 20, units = "in")

arrange association plots

ageeffect_measure = "GAM.smooth.AdjRsq"
ageeffects_plots <- lapply(tracts, plot_age_effect, scalar = "dti_fa", ageeffect_measure = ageeffect_measure, 
                                color_HCPD = color_HCPD, color_HBN = color_HBN, color_PNC = color_PNC, 
                           clipEnds = 5, ylim2 = 0.25, ylim1 = 0, fontsize = 18, legend_position = "none")
names(ageeffects_plots) <- tracts
arrange_association_plots_nolollipop(ageeffects_plots) + bgcolor("white")

#ggsave(paste0(suppfig_dir, "supp_ageeffects_fig2_FA.svg"), arrange_association_plots_nolollipop(ageeffects_plots), height = 10, width = 20, units = "in")
#ggsave(paste0(suppfig_dir, "supp_ageeffects_fig2_FA.png"), arrange_association_plots_nolollipop(ageeffects_plots), height = 10, width = 20, units = "in")

Supplementary Figure 4. Tract-level: age of maturation association with S-A axis.

# reload dti_md files lol
scalar = "dti_md"
HCPD_ageeffects <- read.csv(paste0(output_root, "/", "HCPD", "/GAM/", scalar, "/HCPD_GAM_dev_measures.csv"))
HBN_ageeffects <- read.csv(paste0(output_root, "/", "HBN", "/GAM/", scalar, "/HBN_GAM_dev_measures.csv"))
PNC_ageeffects <- read.csv(paste0(output_root, "/", "PNC", "/GAM/", scalar, "/PNC_GAM_dev_measures.csv"))
# make a list of the different df names: [dataset]_GAM_ageeffects_[scalar]
datasets <- c("HCPD", "HBN", "PNC")
ageeffect_df_names <- c()
for (dataset in datasets) {
  ageeffect_df_names <- append(ageeffect_df_names, paste0(dataset, "_ageeffects"))
}

# format age effect df's
ageeffect_dfs <- lapply(ageeffect_df_names, format_ageeffect)
ageeffect.fdr_dfs <- lapply(ageeffect_dfs, sig_nodes)
## There are 2392/2400 significant nodes
## There are 2350/2400 significant nodes
## There are 2273/2400 significant nodes
names(ageeffect.fdr_dfs) <- ageeffect_df_names
tracts <- setdiff(unique(ageeffect.fdr_dfs$HCPD_ageeffects$tract_label), "Corticospinal") # CST for supplement
# load glasser labels
glasser_labels <- read.csv("/cbica/projects/luo_wm_dev/atlases/glasser/HCP-MMP1_UniqueRegionList.csv")
glasser_labels$regionID <- c(1:360)
glasser_labels$region <- gsub("7Pl", "7PL", glasser_labels$region)   

# load maps for depth = 1.5mm (disregard a,b, and c. Variables are just there to prevent lapply output from being printed lol)
a <- lapply("1.5", load_maps, "PNC") # makes lh_maps_PNC, rh_maps_PNC
b <- lapply("1.5", load_maps, "HCPD") # makes lh_maps_HCPD, rh_maps_HCPD 
c <- lapply("1.5", load_maps, "HBN") # makes lh_maps_HBN, rh_maps_HBN
glasser_SAaxis <- read.csv("/cbica/projects/luo_wm_dev/SAaxis/glasser_SAaxis.csv")
glasser_SAaxis <- glasser_SAaxis %>% select(SA.axis_rank, label)
glasser_SAaxis$regionName <- gsub("_ROI", "", glasser_SAaxis$label)
glasser_SAaxis$regionName <- gsub("^(.)_(.*)$", "\\2_\\1", glasser_SAaxis$region)
threshold=0.3

PNC_deveffects_5_agemat <- ageeffect.fdr_dfs$PNC_ageeffects %>% filter((nodeID < 10 & nodeID > 4) | (nodeID < 95 & nodeID > 89)) %>% mutate(node_position = case_when((nodeID < 10 & nodeID > 4) ~ "end1", (nodeID < 95 & nodeID > 89) ~ "end2")) %>%
  select(tract_label, tractID, nodeID, node_position, hemi, smooth.peak.change, smooth.decrease.offset, smooth.last.change, smooth.slowing.onset) %>%
  group_by(tractID, node_position, hemi) %>% summarise(mean_peak_change = mean(smooth.peak.change, na.rm = T),
                              mean_ageeffect = mean(smooth.decrease.offset, na.rm = T),
                              mean_last_change = mean(smooth.last.change, na.rm = T),
                              mean_dev_slowing = mean(smooth.slowing.onset, na.rm = T))

HCPD_deveffects_5_agemat <- ageeffect.fdr_dfs$HCPD_ageeffects %>% filter((nodeID < 10 & nodeID > 4) | (nodeID < 95 & nodeID > 89)) %>% mutate(node_position = case_when((nodeID < 11 & nodeID > 4) ~ "end1", (nodeID < 95 & nodeID > 89) ~ "end2")) %>%
  select(tract_label, tractID, nodeID, node_position, hemi, smooth.peak.change, smooth.decrease.offset, smooth.last.change, smooth.slowing.onset) %>%
  group_by(tractID, node_position, hemi) %>% summarise(mean_peak_change = mean(smooth.peak.change, na.rm = T),
                              mean_ageeffect = mean(smooth.decrease.offset, na.rm = T), # mean_ageeffect = mean_age_mat
                              mean_last_change = mean(smooth.last.change, na.rm = T),
                              mean_dev_slowing = mean(smooth.slowing.onset, na.rm = T))

HBN_deveffects_5_agemat <- ageeffect.fdr_dfs$HBN_ageeffects %>% filter((nodeID < 10 & nodeID > 4) | (nodeID < 95 & nodeID > 89)) %>% mutate(node_position = case_when((nodeID < 10 & nodeID > 4) ~ "end1", (nodeID < 95 & nodeID > 89) ~ "end2")) %>%
  select(tract_label, tractID, nodeID, node_position, hemi, smooth.peak.change, smooth.decrease.offset, smooth.last.change, smooth.slowing.onset) %>%
  group_by(tractID, node_position, hemi) %>% summarise(mean_peak_change = mean(smooth.peak.change, na.rm = T),
                              mean_ageeffect = mean(smooth.decrease.offset, na.rm = T),
                              mean_last_change = mean(smooth.last.change, na.rm = T),
                              mean_dev_slowing = mean(smooth.slowing.onset, na.rm = T))
       
make_maps("PNC", "5_agemat") # makes [bundle_name]_deveffect_[dataset] for each dataset and bundle. e.g. IFOL_deveffect_PNC
make_maps("HCPD", "5_agemat") 
make_maps("HBN", "5_agemat")
region_all_PNC <- aggregate_age_maps("PNC")  
lh_by_region_PNC <- region_all_PNC[[1]]
rh_by_region_PNC <- region_all_PNC[[2]]

region_all_HCPD <- aggregate_age_maps("HCPD")  
lh_by_region_HCPD <- region_all_HCPD[[1]]
rh_by_region_HCPD <- region_all_HCPD[[2]]

region_all_HBN <- aggregate_age_maps("HBN")  
lh_by_region_HBN <- region_all_HBN[[1]]
rh_by_region_HBN <- region_all_HBN[[2]]
all_endpoints_PNC <- compute_mean_SA("PNC")  
lh_all_endpoints_PNC <- all_endpoints_PNC[[1]]
rh_all_endpoints_PNC <- all_endpoints_PNC[[2]]
all_endpoints_PNC <- all_endpoints_PNC[[3]]

all_endpoints_HCPD <- compute_mean_SA("HCPD")  
lh_all_endpoints_HCPD <- all_endpoints_HCPD[[1]]
rh_all_endpoints_HCPD <- all_endpoints_HCPD[[2]]
all_endpoints_HCPD <- all_endpoints_HCPD[[3]]

all_endpoints_HBN <- compute_mean_SA("HBN")  
lh_all_endpoints_HBN <- all_endpoints_HBN[[1]]
rh_all_endpoints_HBN <- all_endpoints_HBN[[2]]
all_endpoints_HBN <- all_endpoints_HBN[[3]]

combined_df <- bind_rows(all_endpoints_PNC %>% mutate(dataset = "PNC"),
                            all_endpoints_HCPD %>% mutate(dataset = "HCPD"),
                            all_endpoints_HBN %>% mutate(dataset = "HBN"))

# Calculate the averages of mean_SA and age_effect across all dataframes
all_endpoints_avgdatasets <- combined_df %>%
  group_by(bundle_name, end)  %>% summarize(mean_SA = mean(mean_SA, na.rm = TRUE),
                                            age_effect = mean(age_effect, na.rm = TRUE),
                                            .groups = "drop")
PNC_tractlevel_SA_pearsons <- read.csv('/cbica/projects/luo_wm_dev/two_axes_manuscript/output/PNC/suppfigs_spintests/supp_tractlevel_pearson_pspin_allendpoints.csv')
HCPD_tractlevel_SA_pearsons <- read.csv('/cbica/projects/luo_wm_dev/two_axes_manuscript/output/HCPD/suppfigs_spintests/supp_tractlevel_pearson_pspin_allendpoints.csv')
HBN_tractlevel_SA_pearsons <- read.csv('/cbica/projects/luo_wm_dev/two_axes_manuscript/output/HBN/suppfigs_spintests/supp_tractlevel_pearson_pspin_allendpoints.csv')
avgdatasets_tractlevel_SA_pearsons <- read.csv('/cbica/projects/luo_wm_dev/two_axes_manuscript/output/avgdatasets/suppfigs_spintests/supp_tractlevel_pearson_pspin_allendpoints.csv')

SA_age_mat_PNC <- plot_meanSA_by_age_mat("PNC", bquote(paste(italic("r"), " = ", .(round(PNC_tractlevel_SA_pearsons$rho.emp, 2)), ", ", italic("p"[spin]), " < 0.0001"))) # pval=0
SA_age_mat_HCPD <- plot_meanSA_by_age_mat("HCPD", bquote(paste(italic("r"), " = ", .(round(HCPD_tractlevel_SA_pearsons$rho.emp, 2)), ", ", italic("p"[spin]), " < 0.0001"))) # pval=0
SA_age_mat_HBN <- plot_meanSA_by_age_mat("HBN", bquote(paste(italic("r"), " = ", .(round(HBN_tractlevel_SA_pearsons$rho.emp, 2)), ", ", italic("p"[spin]), " = N.S.")))
 
SA_age_mat_avgdatasets <- plot_meanSA_by_age_mat("avgdatasets", bquote(paste(italic("r"), " = ", .(round(avgdatasets_tractlevel_SA_pearsons$rho.emp, 2)), ", ", italic("p"[spin]), " < 0.0001"))) + ggtitle("Average Across Datasets")

SA_age_mat_plot_final <- ggarrange(SA_age_mat_PNC, SA_age_mat_HCPD, SA_age_mat_HBN, SA_age_mat_avgdatasets, ncol = 4, common.legend = T, legend = "bottom") 

x.grob <- textGrob("Mean S-A Axis Rank of Cortical Endpoint", 
                   gp=gpar(col="black", fontsize=24))

y.grob <- textGrob("Age of Maturation (Years)", 
                   gp=gpar(col="black", fontsize=24), rot=90)

grid.arrange(arrangeGrob(SA_age_mat_plot_final, left = y.grob, bottom = x.grob))

#ggsave(paste0(suppfig_dir, "supp_SA_tract_level_pearson_allendpoints.png"), grid.arrange(arrangeGrob(SA_age_mat_plot_final, left = y.grob, bottom = x.grob)), height = 6, width = 22, units = "in")
all_endpoints_PNC_binary <- all_endpoints_PNC
max_value <- max(all_endpoints_PNC_binary$age_effect, na.rm = TRUE)
all_endpoints_PNC_binary$age_effect[all_endpoints_PNC_binary$age_effect == max_value] <- NA
all_endpoints_PNC_binary <- all_endpoints_PNC_binary %>% mutate(maturation_status = ifelse(is.na(age_effect), 0, 1)) 
 
all_endpoints_HCPD_binary <- all_endpoints_HCPD
max_value <- max(all_endpoints_HCPD_binary$age_effect, na.rm = TRUE)
all_endpoints_HCPD_binary$age_effect[all_endpoints_HCPD_binary$age_effect == max_value] <- NA
all_endpoints_HCPD_binary <- all_endpoints_HCPD_binary %>% mutate(maturation_status = ifelse(is.na(age_effect), 0, 1)) # 0 = not matured, 1 = matured
 
all_endpoints_HBN_binary <- all_endpoints_HBN
max_value <- max(all_endpoints_HBN_binary$age_effect, na.rm = TRUE)
all_endpoints_HBN_binary$age_effect[all_endpoints_HBN_binary$age_effect == max_value] <- NA
all_endpoints_HBN_binary <- all_endpoints_HBN_binary %>% mutate(maturation_status = ifelse(is.na(age_effect), 0, 1)) 
 
all_endpoints_avgdatasets_binary <- all_endpoints_avgdatasets
max_value <- max(all_endpoints_avgdatasets_binary$age_effect, na.rm = TRUE)
all_endpoints_avgdatasets_binary$age_effect[all_endpoints_avgdatasets_binary$age_effect == max_value] <- NA
all_endpoints_avgdatasets_binary <- all_endpoints_avgdatasets_binary %>% mutate(maturation_status = ifelse(is.na(age_effect), 0, 1))

PNC_tractlevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/suppfigs_spintests/supp_tractlevel_pearson_ttest_pspin_maturedonly.csv")
HCPD_tractlevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/suppfigs_spintests/supp_tractlevel_pearson_ttest_pspin_maturedonly.csv")
HBN_tractlevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/suppfigs_spintests/supp_tractlevel_pearson_ttest_pspin_maturedonly.csv")
avgdatasets_tractlevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/avgdatasets/tract_profiles/suppfigs_spintests/supp_tractlevel_pearson_ttest_pspin_maturedonly.csv")

SA_age_mat_binary_PNC <- plot_meanSA_by_age_mat("PNC", bquote(paste(italic("r"), " = ", .(round(PNC_tractlevel_SA_ttest_p$rho.emp, 2)), ", ", 
                                                      italic("p"[spin]), " < 0.0001")), binary = TRUE) # pval=0
SA_age_mat_binary_HCPD <- plot_meanSA_by_age_mat("HCPD", bquote(paste(italic("r"), " = ", .(round(HCPD_tractlevel_SA_ttest_p$rho.emp, 2)), ", ", 
                                                      italic("p"[spin]), " = ", .(round(HCPD_tractlevel_SA_ttest_p$p.perm.cor, 3)))), binary = TRUE)
SA_age_mat_binary_HBN <- plot_meanSA_by_age_mat("HBN", bquote(paste(italic("r"), " = ", .(round(HBN_tractlevel_SA_ttest_p$rho.emp, 2)), ", ", 
                                                      italic("p"[spin]), " = ", "N.S.")), binary = TRUE)
SA_age_mat_binary_avgdatasets <- plot_meanSA_by_age_mat("avgdatasets", bquote(paste(italic("r"), " = ", .(round(avgdatasets_tractlevel_SA_ttest_p$rho.emp, 2)), ", ", italic("p"[spin]), " < 0.0001")), binary = TRUE) + ggtitle("Average Across Datasets")

SA_age_mat_binary_plot_final <- ggarrange(SA_age_mat_binary_PNC, SA_age_mat_binary_HCPD, SA_age_mat_binary_HBN, SA_age_mat_binary_avgdatasets, ncol = 4, common.legend = T, legend = "bottom")

x.grob <- textGrob("Mean S-A Axis Rank of Cortical Endpoint", 
                   gp=gpar(col="black", fontsize=24))

y.grob <- textGrob("Age of Maturation (Years)", 
                   gp=gpar(col="black", fontsize=24), rot=90)

grid.arrange(arrangeGrob(SA_age_mat_binary_plot_final, left = y.grob, bottom = x.grob))

#ggsave(paste0(suppfig_dir, "supp_SA_tract_level_pearson_maturedOnly.png"), grid.arrange(arrangeGrob(SA_age_mat_binary_plot_final, left = y.grob, bottom = x.grob)), height = 6, width = 22, units = "in")

Supplementary Figure 5. Parcel-level: age of maturation association with S-A axis.

aggregated_axis_PNC <- merge_SA_parcel("PNC")   # 360 glasser parcels with mean age maturation and SA rank
aggregated_axis_HCPD <- merge_SA_parcel("HCPD") 
aggregated_axis_HBN <- merge_SA_parcel("HBN")  
 
# binarize age of maturation df's to get matured and not-yet-matured regions
# make an average binarized df (average across datasets first then binarize)
lh_PNC_merge <- lh_by_region_PNC %>% rename(PNC_mean_ageeffect = regional_mean_ageeffect)
lh_HCPD_merge <- lh_by_region_HCPD %>% rename(HCPD_mean_ageeffect = regional_mean_ageeffect)
lh_HBN_merge <- lh_by_region_HBN %>% rename(HBN_mean_ageeffect = regional_mean_ageeffect)

merge_temp <- merge(lh_HCPD_merge, lh_HBN_merge, by = "region")
lh_all_datasets_ageMat <- merge(merge_temp, lh_PNC_merge, by = "region")
lh_all_datasets_ageMat <- lh_all_datasets_ageMat %>% # left hemi age of maturation maps
  mutate(regional_mean_ageeffect = rowMeans(select(., HCPD_mean_ageeffect, HBN_mean_ageeffect, PNC_mean_ageeffect), na.rm = TRUE))

rh_PNC_merge <- rh_by_region_PNC %>% rename(PNC_mean_ageeffect = regional_mean_ageeffect)
rh_HCPD_merge <- rh_by_region_HCPD %>% rename(HCPD_mean_ageeffect = regional_mean_ageeffect)
rh_HBN_merge <- rh_by_region_HBN %>% rename(HBN_mean_ageeffect = regional_mean_ageeffect)

merge_temp <- merge(rh_HCPD_merge, rh_HBN_merge, by = "region")
rh_all_datasets_ageMat <- merge(merge_temp, rh_PNC_merge, by = "region")
rh_all_datasets_ageMat <- rh_all_datasets_ageMat %>% # right hemi age of maturation maps
  mutate(regional_mean_ageeffect = rowMeans(select(., HCPD_mean_ageeffect, HBN_mean_ageeffect, PNC_mean_ageeffect), na.rm = TRUE))

lh_all_datasets_ageMat$region <- paste0(lh_all_datasets_ageMat$region, "_L")
rh_all_datasets_ageMat$region <- paste0(rh_all_datasets_ageMat$region, "_R")
combined_df <- rbind(lh_all_datasets_ageMat, rh_all_datasets_ageMat) %>% rename(regionName = region)

binarized_combined_df <- right_join(combined_df, glasser_SAaxis, by = "regionName") # binarize
max_value <- max(binarized_combined_df$regional_mean_ageeffect, na.rm = TRUE)
binarized_combined_df$regional_mean_ageeffect[binarized_combined_df$regional_mean_ageeffect == max_value] <- NA

aggregated_axis_avgdatasets <- binarized_combined_df

# load pvalues
PNC_parcellevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/suppfigs_spintests/supp_parcellevel_pearson_pspin.csv")
HCPD_parcellevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/suppfigs_spintests/supp_parcellevel_pearson_pspin.csv")
HBN_parcellevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/suppfigs_spintests/supp_parcellevel_pearson_pspin.csv")
avgdatasets_parcellevel_SA_ttest_p <- read.csv("/cbica/projects/luo_wm_dev/output/avgdatasets/tract_profiles/suppfigs_spintests/supp_parcellevel_pearson_pspin_avgdatasets.csv")

# remember, this data is using bin = 5 (5 nodes at each end, after clipping 5 nodes already)
annot_text_PNC <- bquote(paste(italic("r"), " = ", .(round(PNC_parcellevel_SA_ttest_p$rho.emp, 2)), ", ", 
                                                      italic("p"[spin]), " < 0.0001"))
annot_text_HCPD <- bquote(paste(italic("r"), " = ", .(round(HCPD_parcellevel_SA_ttest_p$rho.emp, 2)), ", ", 
                                                      italic("p"[spin]), " < 0.0001"))
annot_text_HBN <-  bquote(paste(italic("r"), " = ", .(round(HBN_parcellevel_SA_ttest_p$rho.emp, 2)), ", ", 
                                                      italic("p"[spin]), " = ", .(round(HBN_parcellevel_SA_ttest_p$p.perm, 2))))
annot_text_avgdatasets <-  bquote(paste(italic("r"), " = ", .(round(avgdatasets_parcellevel_SA_ttest_p$rho.emp, 2)), ", ", 
                                                      italic("p"[spin]), " < 0.0001")) 

agg_SA_parcel_PNC <- plot_agg_SA_parcel("PNC", annot_text_PNC, all_parcels = TRUE)
agg_SA_parcel_HCPD <- plot_agg_SA_parcel("HCPD", annot_text_HCPD, all_parcels = TRUE)
agg_SA_parcel_HBN <- plot_agg_SA_parcel("HBN", annot_text_HBN, all_parcels = TRUE)
agg_SA_parcel_avgdatasets <- plot_agg_SA_parcel("avgdatasets", annot_text_avgdatasets, all_parcels = TRUE) + ggtitle("Average Across Datasets")

agg_SA_parcel_plot_final <- ggarrange(agg_SA_parcel_PNC, agg_SA_parcel_HCPD, agg_SA_parcel_HBN, agg_SA_parcel_avgdatasets, ncol = 4)
x.grob <- textGrob("S-A Axis Rank", gp=gpar(col="black", fontsize=20))
y.grob <- textGrob("Age of Maturation (Years)", gp=gpar(col="black", fontsize=20), rot=90)
grid.arrange(arrangeGrob(agg_SA_parcel_plot_final, left = y.grob, bottom = x.grob))

#ggsave(paste0(suppfig_dir, "supp_parcel_level_pearsons_allendpoints.png"), grid.arrange(arrangeGrob(agg_SA_parcel_plot_final, left = y.grob, bottom = x.grob)), height = 6, width = 22, units = "in")

Supplementary Figure 7. Figure 1 with ACT in HBN

color_HCPD = "#3FB8AFFF"
color_HBN = "#FF9E9DFF"
color_PNC = "#cc0468"
scalar = "dti_md"
PNC_ageeffects <- read.csv(paste0(output_root, "/", "PNC", "/GAM/", scalar, "/PNC_GAM_dev_measures.csv"))
HCPD_ageeffects <- read.csv(paste0(output_root, "/", "HCPD", "/GAM/", scalar, "/HCPD_GAM_dev_measures.csv"))
HBN_ageeffects <- read.csv(paste0(output_root, "/", "HBN", "/GAM/", scalar, "/HBN_GAM_dev_measures_withACT.csv"))
# make a list of the different df names: [dataset]_GAM_ageeffects_[scalar]
datasets <- c("PNC", "HCPD", "HBN")
ageeffect_df_names <- c()
for (dataset in datasets) {
  ageeffect_df_names <- append(ageeffect_df_names, paste0(dataset, "_ageeffects"))
}

# format age effect df's
ageeffect_dfs <- lapply(ageeffect_df_names, format_ageeffect)
ageeffect.fdr_dfs <- lapply(ageeffect_dfs, sig_nodes)
## There are 2273/2400 significant nodes
## There are 2392/2400 significant nodes
## There are 2175/2200 significant nodes
names(ageeffect.fdr_dfs) <- ageeffect_df_names 
tracts <- setdiff(unique(ageeffect.fdr_dfs$HCPD_ageeffects$tract_label), "Corticospinal") 
# PNC
arc_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Arcuate_bin5_clip5_1sided.RData")
ifo_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Inferior_Fronto_occipital_bin5_clip5_1sided.RData")
ilf_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Inferior_Longitudinal_bin5_clip5_1sided.RData")
parc_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Posterior_Arcuate_bin5_clip5_1sided.RData")
slf_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Superior_Longitudinal_bin5_clip5_1sided.RData")
vof_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Vertical_Occipital_bin5_clip5_1sided.RData")

cc_af_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Anterior_Frontal_bin5_clip5_1sided.RData")
cc_mot_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Motor_bin5_clip5_1sided.RData")
cc_occ_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Occipital_bin5_clip5_1sided.RData")
cc_orb_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Orbital_bin5_clip5_1sided.RData")
cc_pp_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Posterior_Parietal_bin5_clip5_1sided.RData")
cc_sf_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Superior_Frontal_bin5_clip5_1sided.RData")
cc_sp_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Superior_Parietal_bin5_clip5_1sided.RData")
cc_tm_PNC <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/noUF/Callosum_Temporal_bin5_clip5_1sided.RData")
 
NEST_PNC <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC$pval, cc_mot_PNC$pval, cc_occ_PNC$pval, cc_orb_PNC$pval, cc_pp_PNC$pval, cc_sf_PNC$pval, cc_sp_PNC$pval, cc_tm_PNC$pval, arc_PNC$pval, ifo_PNC$pval, ilf_PNC$pval, parc_PNC$pval, slf_PNC$pval, vof_PNC$pval)), method=c("fdr")), rep("PNC", 14))) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))

# HCPD
arc_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Arcuate_bin5_clip5_1sided.RData")
ifo_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Inferior_Fronto_occipital_bin5_clip5_1sided.RData")
ilf_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Inferior_Longitudinal_bin5_clip5_1sided.RData")
parc_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Posterior_Arcuate_bin5_clip5_1sided.RData")
slf_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Superior_Longitudinal_bin5_clip5_1sided.RData")
vof_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Vertical_Occipital_bin5_clip5_1sided.RData")

cc_af_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Anterior_Frontal_bin5_clip5_1sided.RData")
cc_mot_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Motor_bin5_clip5_1sided.RData")
cc_occ_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Occipital_bin5_clip5_1sided.RData")
cc_orb_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Orbital_bin5_clip5_1sided.RData")
cc_pp_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Posterior_Parietal_bin5_clip5_1sided.RData")
cc_sf_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Superior_Frontal_bin5_clip5_1sided.RData")
cc_sp_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Superior_Parietal_bin5_clip5_1sided.RData")
cc_tm_HCPD <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/NEST/noUF/Callosum_Temporal_bin5_clip5_1sided.RData")
 
NEST_HCPD <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_HCPD$pval, cc_mot_HCPD$pval, cc_occ_HCPD$pval, cc_orb_HCPD$pval, cc_pp_HCPD$pval, cc_sf_HCPD$pval, cc_sp_HCPD$pval, cc_tm_HCPD$pval, arc_HCPD$pval, ifo_HCPD$pval, ilf_HCPD$pval, parc_HCPD$pval, slf_HCPD$pval, vof_HCPD$pval)), method=c("fdr")), rep("HCP-D", 14))) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))

# HBN
arc_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Arcuate_bin5_clip5_1sided.RData")
ifo_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Inferior_Fronto_occipital_bin5_clip5_1sided.RData")
ilf_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Inferior_Longitudinal_bin5_clip5_1sided.RData")
parc_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Posterior_Arcuate_bin5_clip5_1sided.RData")
slf_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Superior_Longitudinal_bin5_clip5_1sided.RData")
vof_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Vertical_Occipital_bin5_clip5_1sided.RData")

cc_af_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Anterior_Frontal_bin5_clip5_1sided.RData")
cc_mot_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Motor_bin5_clip5_1sided.RData")
cc_occ_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Occipital_bin5_clip5_1sided.RData")
cc_orb_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Orbital_bin5_clip5_1sided.RData")
cc_pp_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Posterior_Parietal_bin5_clip5_1sided.RData")
cc_sf_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Superior_Frontal_bin5_clip5_1sided.RData")
cc_sp_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Superior_Parietal_bin5_clip5_1sided.RData")
cc_tm_HBN <- readRDS("/cbica/projects/luo_wm_dev/output/HBN/tract_profiles/NEST/noUF/Callosum_Temporal_bin5_clip5_1sided.RData")
 
NEST_HBN <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_HBN$pval, cc_mot_HBN$pval, cc_occ_HBN$pval, cc_orb_HBN$pval, cc_pp_HBN$pval, cc_sf_HBN$pval, cc_sp_HBN$pval, cc_tm_HBN$pval, arc_HBN$pval,  ifo_HBN$pval, ilf_HBN$pval, parc_HBN$pval, slf_HBN$pval, vof_HBN$pval)), method=c("fdr")), rep("HBN", 14))) %>% setNames(c("tract_label", "pval_fdr_all", "Dataset"))
ageeffect_measure = "GAM.smooth.AdjRsq"
ageeffects_plots <- lapply(tracts, plot_age_effect, scalar = "dti_md", ageeffect_measure = ageeffect_measure, 
                                color_HCPD = color_HCPD, color_HBN = color_HBN, color_PNC = color_PNC, clipEnds = 5, ylim2 = 0.41, ylim1 = 0, fontsize = 20, legend_position = "none")
names(ageeffects_plots) <- tracts
colors <- c("#3FB8AFFF", "#FF9E9DFF", "#cc0468")
names(colors) <- c("HCP-D", "HBN", "PNC")
bin_num_nodes = 5
clipEnds = 5

# prepare the data for making lollipop plots
df_all <- bind_rows(
  ageeffect.fdr_dfs$HCPD_ageeffects %>% mutate(Dataset = "HCP-D"),
  ageeffect.fdr_dfs$HBN_ageeffects  %>% mutate(Dataset = "HBN"),
  ageeffect.fdr_dfs$PNC_ageeffects  %>% mutate(Dataset = "PNC")
)

df_all$Dataset <- factor(df_all$Dataset, levels = c("HCP-D", "PNC", "HBN"))
df_formatted <- format_deep_superficial(df_all, ageeffect_measure, bin_num_nodes, clipEnds)  
df_formatted <- df_formatted[complete.cases(df_formatted), ]
 
# merge with NEST results
NEST <- rbind(NEST_HCPD, NEST_HBN, NEST_PNC)
lollipop_data <- prepare_lollipop_data(df_formatted)   
lollipop_data$Dataset <- factor(lollipop_data$Dataset, levels = c("PNC", "HCP-D", "HBN"))

# plot lollipops
tracts <- setdiff(unique(df_formatted$tract_label), "Corticospinal") # CST for supplement
lollipop_plots <- lapply(tracts, make_lollipop_plot, lollipop_data)
names(lollipop_plots) <- tracts
arrange_callosum_plots(lollipop_plots)

#ggsave(paste0(suppfig_dir, "supp_fig1_withact_lollipop.svg"), arrange_callosum_plots(lollipop_plots), height = 13, width = 20, units = "in")
#ggsave(paste0(suppfig_dir, "supp_fig1_withact_lollipop.png"), arrange_callosum_plots(lollipop_plots), height = 13, width = 20, units = "in")

Supplementary Figure 8. Threshold testing for tract-to-cortex maps (10%, 30%, 50%)

# load my colormap :)
colormap_file <- sprintf("%1$s/code/colormaps/aquamarine.json", proj_root)
colormap <- fromJSON(paste(readLines(colormap_file), collapse=""))
aquamarine <- colorRampPalette(colormap)(15)

# arcuate
lh_arc_0.1 <- plot_threshold_maps(lh_maps_PNC$LeftArcuate, 'left', 0.1)
rh_arc_0.1 <- plot_threshold_maps(rh_maps_PNC$RightArcuate, 'right', 0.1)
arc_0.1 <- ggarrange(lh_arc_0.1, rh_arc_0.1, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")  
arc_0.1 <- annotate_figure(arc_0.1, top = text_grob("Arcuate Fasciculus", color = "black", size = 20))

lh_arc_0.3 <- plot_threshold_maps(lh_maps_PNC$LeftArcuate, 'left', 0.3)
rh_arc_0.3 <- plot_threshold_maps(rh_maps_PNC$RightArcuate, 'right', 0.3)
arc_0.3 <- ggarrange(lh_arc_0.3, rh_arc_0.3, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
arc_0.3 <- annotate_figure(arc_0.3, top = text_grob("Arcuate Fasciculus", color = "black", size = 20))

lh_arc_0.5 <- plot_threshold_maps(lh_maps_PNC$LeftArcuate, 'left', 0.5)
rh_arc_0.5 <- plot_threshold_maps(rh_maps_PNC$RightArcuate, 'right', 0.5)
arc_0.5 <- ggarrange(lh_arc_0.5, rh_arc_0.5, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
arc_0.5 <- annotate_figure(arc_0.5, top = text_grob("Arcuate Fasciculus", color = "black", size = 20))

# cst
lh_cst_0.1 <- plot_threshold_maps(lh_maps_PNC$LeftCorticospinal, 'left', 0.1)
rh_cst_0.1 <- plot_threshold_maps(rh_maps_PNC$RightCorticospinal, 'right', 0.1)
cst_0.1 <- ggarrange(lh_cst_0.1, rh_cst_0.1, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
cst_0.1 <- annotate_figure(cst_0.1, top = text_grob("Corticospinal Tract", color = "black", size = 20))

lh_cst_0.3 <- plot_threshold_maps(lh_maps_PNC$LeftCorticospinal, 'left', 0.3)
rh_cst_0.3 <- plot_threshold_maps(rh_maps_PNC$RightCorticospinal, 'right', 0.3)
cst_0.3 <- ggarrange(lh_cst_0.3, rh_cst_0.3, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
cst_0.3 <- annotate_figure(cst_0.3, top = text_grob("Corticospinal Tract", color = "black", size = 20))

lh_cst_0.5 <- plot_threshold_maps(lh_maps_PNC$LeftCorticospinal, 'left', 0.5)
rh_cst_0.5 <- plot_threshold_maps(rh_maps_PNC$RightCorticospinal, 'right', 0.5)
cst_0.5 <- ggarrange(lh_cst_0.5, rh_cst_0.5, nrow = 2, ncol = 1, common.legend = T, legend = "bottom")
cst_0.5 <- annotate_figure(cst_0.5, top = text_grob("Corticospinal Tract", color = "black", size = 20))

threshold_0.1 <- annotate_figure(ggarrange(arc_0.1 + theme(plot.margin = unit(c(0, 0, 1, 0), "cm")), 
                                           cst_0.1 + theme(plot.margin = unit(c(1, 0, 0, 0), "cm")), 
                                           nrow = 2, ncol = 1), top = text_grob("Threshold = 10%", color = "black", size = 20, face = 'bold'))

threshold_0.3 <- annotate_figure(ggarrange(arc_0.3 + theme(plot.margin = unit(c(0, 0, 1, 0), "cm")), 
                                           cst_0.3 + theme(plot.margin = unit(c(1, 0, 0, 0), "cm")), 
                                           nrow = 2, ncol = 1), top = text_grob("Threshold = 30%", color = "black", size = 20, face = 'bold'))

threshold_0.5 <- annotate_figure(ggarrange(arc_0.5 + theme(plot.margin = unit(c(0, 0, 1, 0), "cm")), 
                                           cst_0.5 + theme(plot.margin = unit(c(1, 0, 0, 0), "cm")), 
                                           nrow = 2, ncol = 1), top = text_grob("Threshold = 50%", color = "black", size = 20, face = 'bold'))

all_threshold_plots <- ggarrange(threshold_0.1, threshold_0.3, threshold_0.5, ncol = 3, nrow = 1)
all_threshold_plots

#ggsave(paste0(suppfig_dir, "supp_thresholds_PNC.png"), all_threshold_plots, height = 10, width = 20, units = "in")

Methods: determining best depth for tract-to-cortex mapping

PNC_lh_GM <- read.csv("/cbica/projects/luo_wm_dev/input/PNC/derivatives/gyral_hops/all_subjects/lh_GM_count.csv")
PNC_rh_GM <- read.csv("/cbica/projects/luo_wm_dev/input/PNC/derivatives/gyral_hops/all_subjects/rh_GM_count.csv")
PNC_lh_WM <- read.csv("/cbica/projects/luo_wm_dev/input/PNC/derivatives/gyral_hops/all_subjects/lh_WM_count.csv")
PNC_rh_WM <- read.csv("/cbica/projects/luo_wm_dev/input/PNC/derivatives/gyral_hops/all_subjects/rh_WM_count.csv")

# 1.5 mm is best 
kbl(PNC_lh_GM, align = "lccrr", caption = "PNC Left Hemi: best depths for GM") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")
PNC Left Hemi: best depths for GM
Depth count
depth_0 0
depth_0p5 0
depth_1 375
depth_1p5 726
depth_2 0 
kbl(PNC_rh_GM, align = "lccrr", caption = "PNC Right Hemi: best depths for GM") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")
PNC Right Hemi: best depths for GM
Depth count
depth_0 0
depth_0p5 0
depth_1 385
depth_1p5 716
depth_2 0 
kbl(PNC_lh_WM, align = "lccrr", caption = "PNC Left Hemi: best depths for WM") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")
PNC Left Hemi: best depths for WM
Depth count
depth_0 0
depth_0p5 0
depth_1 272
depth_1p5 829
depth_2 0 
kbl(PNC_rh_WM, align = "lccrr", caption = "PNC Right Hemi: best depths for WM") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")
PNC Right Hemi: best depths for WM
Depth count
depth_0 0
depth_0p5 0
depth_1 273
depth_1p5 828
depth_2 0

Supplementary Table 2: Methods: Deep-Peripheral NEST with different bin sizes (bin size = 3, 5, 7, 10)

tracts <- c("Callosum Anterior Frontal", "Callosum Motor", "Callosum Occipital", "Callosump Orbital", "Callosum Posterior Parietal", "Callosum Superior Frontal", "Callosum Superior Parietal", "Callosum Temporal", "Arcuate", "Inferior Fronto-occipital", "Inferior Longitudinal", "Posterior Arcuate", "Superior Longitudinal", "Uncinate", "Vertical Occipital")

# PNC
# bin size = 3
arc_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Arcuate_bin3_clip5_1sided.RData")
ifo_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Fronto_occipital_bin3_clip5_1sided.RData")
ilf_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Longitudinal_bin3_clip5_1sided.RData")
parc_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Posterior_Arcuate_bin3_clip5_1sided.RData")
slf_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Superior_Longitudinal_bin3_clip5_1sided.RData")
unc_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Uncinate_bin3_clip5_1sided.RData")
vof_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Vertical_Occipital_bin3_clip5_1sided.RData")

cc_af_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Anterior_Frontal_bin3_clip5_1sided.RData")
cc_mot_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Motor_bin3_clip5_1sided.RData")
cc_occ_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Occipital_bin3_clip5_1sided.RData")
cc_orb_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Orbital_bin3_clip5_1sided.RData")
cc_pp_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Posterior_Parietal_bin3_clip5_1sided.RData")
cc_sf_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Frontal_bin3_clip5_1sided.RData")
cc_sp_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Parietal_bin3_clip5_1sided.RData")
cc_tm_PNC_bin3 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Temporal_bin3_clip5_1sided.RData")
 
NEST_PNC_bin3 <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC_bin3$pval, cc_mot_PNC_bin3$pval, cc_occ_PNC_bin3$pval, cc_orb_PNC_bin3$pval, cc_pp_PNC_bin3$pval, cc_sf_PNC_bin3$pval, cc_sp_PNC_bin3$pval, cc_tm_PNC_bin3$pval, arc_PNC_bin3$pval, ifo_PNC_bin3$pval, ilf_PNC_bin3$pval, parc_PNC_bin3$pval, slf_PNC_bin3$pval, unc_PNC_bin3$pval, vof_PNC_bin3$pval)), method=c("fdr")), rep("PNC", 15))) %>% setNames(c("tract_label", "bin_3_pvalue", "Dataset"))
 
# bin size = 5
arc_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Arcuate_bin5_clip5_1sided.RData")
ifo_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Fronto_occipital_bin5_clip5_1sided.RData")
ilf_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Longitudinal_bin5_clip5_1sided.RData")
parc_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Posterior_Arcuate_bin5_clip5_1sided.RData")
slf_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Superior_Longitudinal_bin5_clip5_1sided.RData")
unc_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Uncinate_bin5_clip5_1sided.RData")
vof_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Vertical_Occipital_bin5_clip5_1sided.RData")

cc_af_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Anterior_Frontal_bin5_clip5_1sided.RData")
cc_mot_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Motor_bin5_clip5_1sided.RData")
cc_occ_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Occipital_bin5_clip5_1sided.RData")
cc_orb_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Orbital_bin5_clip5_1sided.RData")
cc_pp_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Posterior_Parietal_bin5_clip5_1sided.RData")
cc_sf_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Frontal_bin5_clip5_1sided.RData")
cc_sp_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Parietal_bin5_clip5_1sided.RData")
cc_tm_PNC_bin5 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Temporal_bin5_clip5_1sided.RData")
 
NEST_PNC_bin5 <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC_bin5$pval, cc_mot_PNC_bin5$pval, cc_occ_PNC_bin5$pval, cc_orb_PNC_bin5$pval, cc_pp_PNC_bin5$pval, cc_sf_PNC_bin5$pval, cc_sp_PNC_bin5$pval, cc_tm_PNC_bin5$pval, arc_PNC_bin5$pval, ifo_PNC_bin5$pval, ilf_PNC_bin5$pval, parc_PNC_bin5$pval, slf_PNC_bin5$pval, unc_PNC_bin5$pval, vof_PNC_bin5$pval)), method=c("fdr")), rep("PNC", 15))) %>% setNames(c("tract_label", "bin_5_pvalue", "Dataset"))
 
# bin size = 7
arc_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Arcuate_bin7_clip5_1sided.RData")
ifo_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Fronto_occipital_bin7_clip5_1sided.RData")
ilf_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Longitudinal_bin7_clip5_1sided.RData")
parc_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Posterior_Arcuate_bin7_clip5_1sided.RData")
slf_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Superior_Longitudinal_bin7_clip5_1sided.RData")
unc_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Uncinate_bin7_clip5_1sided.RData")
vof_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Vertical_Occipital_bin7_clip5_1sided.RData")

cc_af_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Anterior_Frontal_bin7_clip5_1sided.RData")
cc_mot_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Motor_bin7_clip5_1sided.RData")
cc_occ_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Occipital_bin7_clip5_1sided.RData")
cc_orb_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Orbital_bin7_clip5_1sided.RData")
cc_pp_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Posterior_Parietal_bin7_clip5_1sided.RData")
cc_sf_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Frontal_bin7_clip5_1sided.RData")
cc_sp_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Parietal_bin7_clip5_1sided.RData")
cc_tm_PNC_bin7 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Temporal_bin7_clip5_1sided.RData")
 
NEST_PNC_bin7 <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC_bin7$pval, cc_mot_PNC_bin7$pval, cc_occ_PNC_bin7$pval, cc_orb_PNC_bin7$pval, cc_pp_PNC_bin7$pval, cc_sf_PNC_bin7$pval, cc_sp_PNC_bin7$pval, cc_tm_PNC_bin7$pval, arc_PNC_bin7$pval, ifo_PNC_bin7$pval, ilf_PNC_bin7$pval, parc_PNC_bin7$pval, slf_PNC_bin7$pval, unc_PNC_bin7$pval, vof_PNC_bin7$pval)), method=c("fdr")), rep("PNC", 15))) %>% setNames(c("tract_label", "bin_7_pvalue", "Dataset"))
 
# bin size = 10
arc_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Arcuate_bin10_clip5_1sided.RData")
ifo_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Fronto_occipital_bin10_clip5_1sided.RData")
ilf_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Inferior_Longitudinal_bin10_clip5_1sided.RData")
parc_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Posterior_Arcuate_bin10_clip5_1sided.RData")
slf_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Superior_Longitudinal_bin10_clip5_1sided.RData")
unc_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Uncinate_bin10_clip5_1sided.RData")
vof_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Vertical_Occipital_bin10_clip5_1sided.RData")

cc_af_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Anterior_Frontal_bin10_clip5_1sided.RData")
cc_mot_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Motor_bin10_clip5_1sided.RData")
cc_occ_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Occipital_bin10_clip5_1sided.RData")
cc_orb_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Orbital_bin10_clip5_1sided.RData")
cc_pp_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Posterior_Parietal_bin10_clip5_1sided.RData")
cc_sf_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Frontal_bin10_clip5_1sided.RData")
cc_sp_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Superior_Parietal_bin10_clip5_1sided.RData")
cc_tm_PNC_bin10 <- readRDS("/cbica/projects/luo_wm_dev/output/PNC/tract_profiles/NEST/Callosum_Temporal_bin10_clip5_1sided.RData")
 
NEST_PNC_bin10 <- data.frame(cbind(tracts, p.adjust(unlist(c(cc_af_PNC_bin10$pval, cc_mot_PNC_bin10$pval, cc_occ_PNC_bin10$pval, cc_orb_PNC_bin10$pval, cc_pp_PNC_bin10$pval, cc_sf_PNC_bin10$pval, cc_sp_PNC_bin10$pval, cc_tm_PNC_bin10$pval, arc_PNC_bin10$pval, ifo_PNC_bin10$pval, ilf_PNC_bin10$pval, parc_PNC_bin10$pval, slf_PNC_bin10$pval, unc_PNC_bin10$pval, vof_PNC_bin10$pval)), method=c("fdr")), rep("PNC", 15))) %>% setNames(c("tract_label", "bin_10_pvalue", "Dataset"))
PNC_bin_compare <- left_join(NEST_PNC_bin3, NEST_PNC_bin5) %>% left_join(NEST_PNC_bin7) %>% left_join(NEST_PNC_bin10) %>% select(-Dataset)
PNC_bin_compare$bin_3_pvalue <- as.numeric(PNC_bin_compare$bin_3_pvalue)
PNC_bin_compare$bin_5_pvalue <- as.numeric(PNC_bin_compare$bin_5_pvalue)
PNC_bin_compare$bin_7_pvalue <- as.numeric(PNC_bin_compare$bin_7_pvalue)
PNC_bin_compare$bin_10_pvalue <- as.numeric(PNC_bin_compare$bin_10_pvalue)
# binning differences don’t impact significance status, though specific p-value changes a little
PNC_bin_table <- PNC_bin_compare %>%
  mutate(across(where(is.numeric), ~ formatC(.x, digits = 2, format = "fg", drop0trailing = TRUE))) %>%
  rename(Tract = tract_label,
         "Bin Size 3" = bin_3_pvalue,
         "Bin Size 5" = bin_5_pvalue,
         "Bin Size 7" = bin_7_pvalue,
         "Bin Size 10" = bin_10_pvalue)

PNC_bin_table_final <- PNC_bin_table %>%
  kbl(format = "html", align = "lcccc", digits = 2) %>%
  add_header_above(c(" " = 1, "p-value" = 4), bold = TRUE) %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(0:nrow(PNC_bin_table), extra_css = "line-height: 0.5;")  

PNC_bin_table_final
p-value
Tract Bin Size 3 Bin Size 5 Bin Size 7 Bin Size 10
Callosum Anterior Frontal 0.00012 0.00012 0.00012 0.00014
Callosum Motor 0.00012 0.00012 0.00012 0.00014
Callosum Occipital 0.00012 0.00012 0.00012 0.032
Callosump Orbital 0.00012 0.00012 0.00012 0.00014
Callosum Posterior Parietal 0.00012 0.00012 0.00012 0.00014
Callosum Superior Frontal 0.00012 0.00012 0.00012 0.00014
Callosum Superior Parietal 0.00012 0.00012 0.00012 0.00014
Callosum Temporal 0.00012 0.00012 0.00012 0.00014
Arcuate 0.00012 0.00012 0.00012 0.00014
Inferior Fronto-occipital 0.00012 0.00012 0.0093 0.0075
Inferior Longitudinal 0.36 0.22 0.26 0.24
Posterior Arcuate 0.00012 0.00012 0.00012 0.00014
Superior Longitudinal 0.00012 0.00012 0.00012 0.00014
Uncinate 0.17 0.15 0.12 0.1
Vertical Occipital 0.00012 0.00012 0.00012 0.00014 
#save_kable(PNC_bin_table_final, paste0(suppfig_dir, "table2_binning.html"))
#webshot(paste0(suppfig_dir, "table2_binning.html"), paste0(suppfig_dir, "table2_binning.png"))
#webshot(paste0(suppfig_dir, "table2_binning.html"), paste0(suppfig_dir, "table2_binning.pdf"))